clear
close all
clc

ln = @(x)(log(x));

%% Parametros de ajuste de la simulacion
M = 2;  % cantidad de simbolos
k = log2(M); % bits por simbolo
SNRbdB = 2:14; % SNR por bit en dB
nsimbs = 2e6; % cantidad de simbolos que para cada valor de SNRbdB en el Montercarlo


%% Conversiones e inicializaciones varias
SNRbveces = 10.^(SNRbdB/10); % SNR por bit en veces
SNRdB = SNRbdB + 10 * log10(k);  % SNR en dB
SNR = 10.^(SNRdB/10);    % SNR en veces
d = 2;
Eavg = (M^2-1)/12 * d^2;

Pes = zeros(size(SNRbdB));  % probabilidad de error por simbolo simulada
Peb = zeros(size(SNRbdB));  % probabilidad de error de bit simulada

%No importa el valor, pues en la expresion de la Pe se anula su efecto
No = 1;

p = 1/4;
A = 1;
%Threshhold
Th = 1/(4*A) * No * ln ((1-p)/p);

%% P(error) teorica:
% Probabilidad de error de bit teorica de una 4-PAM
% Uso un puntero a una funcion cuyo parametro es SNR (en veces)
Pet = @(SNR) p * (1 - qfunc((Th - A)/(No/2))) + (1 - p) * qfunc((Th + A)/(No/2));

%Pet = @(SNR) 2 * (1-1/M) * qfunc(sqrt(3 * SNR/(M^2-1)));

%% Simulacion Montecarlo

Th_vector = ones(nsimbs,1) * Th; 

for n0 = 1:length(SNRbdB)
    
    %% Primero trabajo con simbolos y despues analizo por bit
    
    sigma = sqrt(Eavg/SNR(n0)); % desvio de ruido
    %simbtx = randi(M,nsimbs,1)-1;   % simbolos a transmitir (en 1:M);
    
    simbtx = binornd(1, p,nsimbs,1);
   
    % Sumo el ruido (OJO sumo el ruido como si la constelacion tuviera parte compleja)
    simbmod = simbtx + (sigma) * (randn(nsimbs,1)+1i * randn(nsimbs,1));
         
    % Probabilidad de error de simbolo simulada
    Pes(n0) = nnz(simbtx>=Th_vector)/nsimbs;
    
end


%% Graficos
% Probabilidad de error de simbolo
figure(1)
semilogy(SNRbdB,Pet(SNR),'-dr','linewidth',1.5);
hold on
semilogy(SNRbdB,Pes,'-o','linewidth',1.5);
grid on
box on
xlabel('SNR por bit (dB)');
ylabel('P_e')
legend('Curva teorica','Montecarlo')

